---
title: "Fast expansion of a polynomial with Julia"
author: "Stéphane Laurent"
date: '2022-07-17'
tags: julia, povray, maths
output:
  md_document:
    variant: markdown
    preserve_yaml: true
  html_document:
    highlight: kate
    keep_md: no
highlighter: pandoc-solarized
---

In the [previous
post](https://laustep.github.io/stlahblog/posts/caracas02.html), I
showed how to expand a polynomial in R when the coefficients of the
polynomial contain some literal values, with the help of the **caracas**
package (a wrapper of **SymPy**). As I wrote, the method takes 20
minutes for the example I'm interested in.

Then I tried to do it with Julia. My attempt has been successful, and as
usual Julia astonished me: it took about one minute only!

Here is the code, for my real example.

``` {.julia .numberLines}
using AbstractAlgebra
using Printf # used to write the result

# define the symbolic constants
SS, (a,b,c,d,sqrt3,w0) = PolynomialRing(QQ, ["a","b","c","d","sqrt3","w0"])
# define the symbolic variables
TT, (x,y,z,w) = PolynomialRing(SS, ["x","y","z","w"])

# polynomial expression
P = ((x^2+y^2+z^2+w^2+145//3)^2-4*(9*z^2+16*w^2))^2*((x^2+y^2+z^2+w^2+145//3)^2+296*(x^2+y^2)-4*(9*z^2+16*w^2))-16*(x^2+y^2)*(x^2+y^2+z^2+w^2+145//3)^2*(37*(x^2+y^2+z^2+w^2+145//3)^2-1369*(x^2+y^2)-7*(225*z^2+448*w^2))-sqrt3*16//9*(x^3-3*x*y^2)*(110*(x^2+y^2+z^2+w^2+145//3)^3-148*(x^2+y^2+z^2+w^2+145//3)*(110*x^2+110*y^2-297*z^2+480*w^2))-64*(x^2+y^2)*(3*(729*z^4+4096*w^4)+168*(x^2+y^2)*(15*z^2-22*w^2))+64*(12100//27*(x^3-3*x*y^2)^2-7056*(3*x^2*y-y^3)^2)-592240896*z^2*w^2

# substitutions to be done
X = a*x - b*y - c*z - d*w0
Y = a*y + b*x + c*w0 - d*z
Z = a*z - b*w0 + c*x + d*y
W = a*w0 + b*z - c*y + d*x

# perform the substitutions
Q = P(X, Y, Z, W)

# vector of powers (exponents)
pows = map(degrees, terms(Q))
# vector of coefficients as strings
coeffs_str = map(string, coefficients(Q))

# generate the POV-Ray code
# three replacements:
# fractions: 145//3 => 145/3
# powers: w0^2 => pow(w0,2), a^3 => pow(a,3), b^4 => pow(b,4)
function g(degrees_coeff)
  d, c = degrees_coeff
  c = replace(c, r"//" => s"/")
  c = replace(c, r"w0\^(\d+)" => s"pow(w0,\g<1>)")
  c = replace(c, r"([ab])\^(\d+)" => s"pow(\g<1>,\g<2>)")
  return @sprintf "xyz(%d,%d,%d): %s,\n" d[1] d[2] d[3] c
end
povray = map(g, zip(pows, coeffs_str))

# write the POV-Ray code to a file
io = open("povray.txt", "w")
for trm in povray
    write(io, trm)
end
close(io)
```

Not only it is fast, but it is also more easy to read than the R/Python
code.

Here is how the generated file **povray.txt** looks like (it takes more
than one megabyte):

![](./figures/povray_polynomial.png)